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ABSTRACT 

The past five years has seen a surge in research and innovative ideas for the imaging of 
extrasolar planets, particular terrestrial ones. We expect that within the next decade a space ob- 
servatory will be launched with the objective of imaging earthlike planets. Because of the limited 
lifetime of such a mission and the large number of potential targets, integration time is a critical 
parameter. In fact, integration time is the primary metric in evaluating various design approaches 
for the high contrast imaging system. In this paper we present a new approach to determining 
the existence of a planet in an observed system using Bayesian hypothesis testing. Rather than 
perform photometry, or rely on vision to determine the existence of a planet, this approach eval- 
uates the image plane data statistically under certain assumptions about the prior probability 
distributions. We show that extremely high confidence can be achieved in substantially shorter 
integration times than conventional photometric methods. 



1. Introduction 

The discovery of more than 100 extrasolar Jupiter-sized planets in just the last decade has generated 
enormous interest, both among astronomers and the public, in the problem of discovering and characterizing 
Earthlikc planets. NASA is already planning its next large space-based observatory, the Terrestrial Planet 
Finder (TPF), with a planned launch date toward the end of the next decade. TPF's primary objective will 
be to discover Earthlike planets and characterize them for indications of life. 

The technical challenges for TPF are great. Foremost among them is the problem of high-contrast 
imaging. In order to discover as many planets as possible, it is necessary to design an imaging systems 
that achieves very high contrast between the parent star and the planet as close as possible to the star. An 
earlier study by Brown et al. (2002) indicates that a D = 4m class visible-light (i.e. 400nm < A < 650nm) 
instrument ought to be able to discover about 50 extrasolar Earth-like planets if it can provide contrast of 
10~ 10 at an angular separation of 3X/D and that a 4 x 10m class telescope ought to be able to discover 
about 150 such planets if it can provide the same contrast at a separation of 4A/ D. 



1 I. Braems received a fellowship from the Institut National de Recherche en Informatique et Automatique (INRIA) 
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There are myriad of methods being proposed to achieve the high contrast necessary for planet detection. 
Our interest has been in shaped pupil coronagraphy. Many of our optimal approaches are described in Kasdin 
ct al. (2003); Vanderbei et al. (2003a, b). Of particular importance is the discussion in Kasdin et al. (2003) 
on performance metrics used to compare and optimize various coronagraphic approaches. The key factor 
here is integration time. All planet finding missions must minimize the integration time in order to provide 
the most possible observations and the least sensitivity to instabilities in the telescope and spacecraft. To 
that end, quantitative measures of integration time, derived for specific image processing techniques, are 
essential. 

In Kasdin ct al. (2003) we discussed a number of different approaches for planet detection, each with 
a corresponding integration time to reach a specified signal-to-noise ratio (S/N). The objective was to find 
quantitative metrics that could be used to compare different coronagraph approaches. The simplest detection 
integration time threshold, t±, was proposed such that for r > ti the decision of whether a planet exists can 
be taken by a human eye. This threshold is based on the expected signal to noise ratio S/N 

where it is assumed that the planet intensity I p is known, the background is uniform and its intensity lb 
is also know and AP is the area under the PSF P in the integration area AS around the planet location 
to be denoted £* , where £* refers to the (continuous) location in the image plane of the central peak of 
the planet's PSF (note that C* can be considered a member of M 2 , thus representing the two-dimensional 
position of the planet on the image plane). The usual assumption is that the human eye can detect a planet 
when the S/N ratio exceeds 5, setting the typical detection criterion. Note the dependence on PSF shape; 
as the coronagraphs we are considering operate by modifying the shape of the PSF, there is some sensitivity 
of the integration time to the specific coronagraph design. Brown and Burrows (1990) call this ratio of PSF 
area to pixel area the "sharpness" . 

Unfortunately, this photometric approach to detection is limited, as it relies on an empirical criterion 
that does not provide any information about the confidence level to be associated with the collected images. 
In Kasdin ct al. (2003), other integration times were computed but they were based on the assumption 
that photometry was performed on the image in an effort to estimate the planet (and perhaps background) 
irradiance. The goal of this paper is to establish a simple detection criterion to be used before any photom- 
etry analysis in order to allow a smaller integration time, while providing quantitative information on the 
confidence associated with a planet decision. We also plan to study modifications of the coronagraph design 
optimization using these results. 

The basic principal is to use probabilistic hypothesis testing to confirm or deny the hypothesis that a 
planet is located at a given pixel. Among several hypothesis testing techniques, Bayesian approaches seem 
promising. Their appeal stems, in large part, from the fact that these tests are conditioned on the collected 
data, no matter the amount, while frequentist approaches require interpolation on fictitious data (as multiple 
data sets for determining statistics will not be available). Moreover, the mathematical framework of Bayesian 
methods includes marginalization, i.e., an easy way to take into account uncertain nuisance parameters (see 
Berger (1985)). For these reasons and others, Bayesian techniques have resulted in a significant breakthrough 
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in the last decade (see for example Loredo (1992), Gregory and Loredo (1992), Defay et al. (2001) and Aigrain 
and Favata (2002)). 

For the sake of simplicity and demonstration, we only consider ID images in this work (that is, the 
dependence of the PSF along only one axis of the image plane). While the real image is a continuous 
function of the position (: T (C) = h + I P P (C — C*)> the CCD collects a photon noise corrupted discretized 
image Z = (z n ) n _ 1 M where z n stands for the intensity of the n — th pixel. We use n* to indicate the index 
of the pixel containing (*: (* e [n'Aa, (n* + l)Aa]. Due to the low photon count at these short integration 
times, each z n is modeled as a Poisson random variable whose mean is 

X n (I b , I p , C) = r (^I b Aa + I p J^ +1 P (C - O <Kj = r (I b Aa + I p AP n (O) , Vn - 1, M. (2) 

where Aa is the size of a pixel and (£„, Cn+i) refer to the location of the edges of pixel n. In this preliminary 
study, we will assume that lb is uniform, we will neglect the speckle noise due to imperfect optics, and we 
will consider only monochromatic signals. Clearly these are rather restrictive assumptions, but they allow 
us to focus attention on the technique and compare our results to the previous one. We hope during the 
next few years to expand this effort and use the same Bayesian ideas on the more realistic polychromatic 
system with speckle. For instance, a planet might be separated from speckle by incorporating into the model 
a few discrete frequencies (most likely via a dichroic) or via multiple measurements at different observatory 
orientations. 



2. Detection criterion 

In most cases (except Gregory and Loredo (1992)) Bayesian approaches were used only to solve the 
inverse problem of denoising the collected image. This would correspond to estimating I p (i.e., it corresponds 
to the photometric approaches proposed in Kasdin et al. (2003)) or to estimating (* (the localization problem) 
from Eq. 2. Both of these approaches by necessity already assume that a planet is present and typically take 
more time. A hypothesis testing approach is a direct approach that will more appropriately be performed 
before any further processing. 



2.1. Odds ratio 

As the objective is not only to detect a potential planet but also its discrete location n* , we propose a 
test to be performed on each pixel. To determine whether a potential planet is present in the n — th pixel, 
we compare the probabilities of the following two alternate hypotheses: 

Hq : there is no planet in the n — th pixel, 
H" : there is a planet in the n — th pixel. 
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Associate a model M, to each hypothesis Hi: 

M " : \ n = \ n (I b ,I p ,(*), C* € [nAa, (n + l)Aa], J p = 
M'\ : X n = X n (I b ,I p ,C), C € [nAa,(n + l)Aa], J p ^ 

We now compare the posterior probabilities of the two models conditioned on the collected data, and 
we shall decide that there is a planet in pixel n if the odds ratio, O™ , favors iff over iff: 

= p (Mi\Z) p(M?)p(Z\M?) 
10 p(M »|Z) p(M «)p(Z|M «) 1 ; 

where we have used Bayes' rule to write the odds ratio as the product of the ratio of prior probabilities ^pjrrj 

and the ratio of likelihoods ^t^tjjX- This dependence on the prior probability of a terrestrial planet in a 

given system (often referred to as r/Earth) is one of the most attractive features of the Bayesian approach, as 
information from prior missions such as Kepler can be explicitly included in the detection algorithm and the 
sensitivity of the detection to these assumptions can be quantitatively assessed. 

One possible decision criterion would be that O™ > 1 favors if™, thus indicating that there is a planet 
whose true location is inside the n — th pixel. Alternatively, the decision can be made based on a different 
value of the ratio derived from particular weighting criteria. We call this a "Loss-based decision" and describe 
it next. 



2.2. Loss-based decision 

As mentioned above, a decision criterion based on an odds ratio of unity is not the only possible one. 
In particular, an alternative criteria uses a weighted decision based upon our sensitivity to different types of 
errors (for instance, we may be more willing to tolerate missed detections than false alarms). To quantify 
this we define two possible actions: 

a : "we accept if " a\ : "we reject ifo" . 

The Receiving Operator Curve (ROC) is used to determine the best possible decision threshold and quantify 
the relative importance of the two types of errors. The ROC is a plot of the evolution of the missed detection 
probability (Type II error), p(a |ff™), versus 1 — p(a\\HQ ), where p{cli\Hq ) is the false alarm rate (Type 
I error). In other words, we place a point on the curve corresponding to the two probabilities for every 
possible threshold value, from to infinity. We then typically select the threshold corresponding to the 
"knee" of the curve, that is, the combination of maximum probability of success with minimum likelihood 
of missed detection. The area under the ROC is defined as the "power" of the test and for the results below 
is parameterized by the integration time. A perfect test is one with zero probability of missed detection and 
zero probability of false alarm, resulting in a power value (area under the ROC) of one. The values of the 
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probabilities for the curve are found through Monte Carlo simulations. 2 

2.3. Likelihoods and priors 

To compute O™ as defined in Eq. 3, we need to assign the priors p(M") , i = 1,2 and compute the 
likelihoods p (Z|M™) , i = 1, 2. In the absence of any other information, we assume here that the two models 
are a priori equiprobable at each pixel: p (Mq) = \, p {M™) = 5, Vn = 1, . . . , M. 

In general, the probability of a given data set depends on various physical parameters of the system. 
For instance, the mean photon rate in the Poisson distribution depends on the irradiance and the planet 
location in the pixel. In order to compute the likelihood, we must marginalize the probability to remove this 
parameter dependence. The marginalization law allows us to compute the global likelihood of the model 
M™(ip) over its parameter vector denoted by tp: 

P (z\m?) = J d^ P (m?)p(z\^Mn- (4) 

Modeling each pixel intensity Zj as an independent Poisson random variable, we obtain 

M M . , . , Zj 

P m,M?) = I]p(*ihW) = n , ' exp-(A(^M?) 3 -) (5) 

3=1 j=i Zj ' 

which we multiply by the prior probabilities and then integrate over ip to determine the likelihoods as in 
Eq. 4. 

To determine the prior probabilities, we assume that all the parameters are independent, resulting in: 

p (I p ,I b , C\M?) = p (I p \M?)p(I b \M?)p(C\M?) , 

where here ip consists of the planet and background irradiances and the planet location. We thus need to 
assign prior probabilities to each parameter individually. The choice of specific densities for the priors is 
a delicate subject as there is not necessarily a sound theoretical basis for the choice and particular priors 
may effect the performance of the estimator. Our approach here is to make a sensible choice for the priors 
that also makes the mathematical analysis manageable. In future work we will examine the effect on the 
probability of success when the assumed prior densities are actually in error. 

The continuous planet location £* is assumed to be uniform over a pixel containing a planet: 

p((*\M?) = -|- when C* e [nAa, (n+l)Aa], 
= otherwise. 



2 It was found empirically that the optimal threshold value strongly depends upon the shape of the PSF, with a unity 
threshold best only for a square PSF. 
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As the intensity I p of the planet is strictly positive, it has been modeled by a Gamma law T (a, [3) 
of mean E[I P ] = a(3 and variance E[I%] — a/3 2 . A second advantage of this model is that the gamma 
distribution is the conjugate prior of a Poisson law (Berger (1985)), thus simplifying future computations: 

P (/ p |M 1 ") = /;- lex p(-^)- ( fi ) 

We will either assume that the intensity of the background is known or that it also follows a Gamma 
law, T (a, b). 



3. Results 

For simplicity all the results provided in this section have been performed for a 1-D point spread function, 
that is, a cross section of the two-dimensional image. This function for optimal shaped pupil coronagraphs is 
the transform of the prolate spheroidal wavefunction (see Kasdin et al. (2003) or Vanderbei et al. (2003a)). 
In future work we will repeat and compare to a number of different pupil designs (as well as other types of 
coronagraphs) and extend to two-dimensions. All results are based on Monte-Carlo simulations. 

To account for the size of the PSF, which extends over more than one pixel, we define A as the width 
of the significant part of the PSF (normally taken equal to the FWHM), and set 

such that 2K + 1 is the maximum number of pixels that can be covered by the truncated PSF and AS = 
(2K + l)Aa. For a pixel nj such that \\rij - n\\ > K, A(V>, M^- = A(^,Mf)j = I b Aar . The ratio of 
likelihoods as defined in Eq. 5 only involves the pixels covered by A: 



I? 

ram(n+K,M) 



p(zw,mz) njiiPtehw) 



n 



mi „(„ + *,M) /x{ ^ m 



n 



j — max ( n —K,l) 

where M is the total number of pixels. 



A(V,M ") 



exp-(A(^,M 1 "),-A(^M ") i ) 



In the following simulations we have selected the pixel size such that K = 2, and for the sake of simplicity, 
we will consider n such that max(n — K, 1) = n — K, min(n + K, M) = n + K (no edge). 
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3.1. Assume I, p and I(, are known 

For this study, we use Monte-Carlo simulations to determine the effectiveness of the hypothesis testing 
criterion and the resulting probabilities of type I and II errors for a given normalized integration time, r. 
This can then be compared to the decision criterion described in Kasdin et al. (2003) and based on Eq. 1 for 
<i = t. In this section, the simulations assume I p and lb are both known. In reality, of course, not only is I p 
not necessarily known, it varies with pixel location as the planet moves further away from the star (assuming 
a reflected visible light system). This assumption is relaxed in later sections. 

Under these assumptions, the odds ratio in favor of H[ l is (see Appendix): 

r Aa 3=n+K . . Zj 

O? = —exp(-Tl p AP) [] 1+ A^ AP ^ <W ' (7) 

J ° j=n-K ^ ' 

where Q — and APj(9) is the area of the PSF covered by one of the five pixels (indexed by j) when 
the planet appears at location 8 in pixel n. Figure 1 depicts the results obtained for two simulated sets 
of images. For i e [1,500], the data are simulated without any planet (Q = 0), while for i 6 [501, 1000], 
Q = 1, the planet is located in pixel n* = 11. For all simulations, the normalized integration time was set 
to t = 1, corresponding to a signal-to-noise ratio S/N computed as in Eq. 1 of 3.52. While it is reasonable 
to expect some detections using the conventional approach with this S/N ratio, using the usual criteria of 
S/N = 5 described above, no planets would be detected. Figure 2 depicts the histogram of logO™ for the 
two simulated sets. For a threshold equal to one, the false alarm rate is 0.1%, while the missed detection rate 
reaches 7.2%. From the ROC curve (not plotted here for brevity), we determined that the test has maximal 
power if the threshold is set to 0.158 rather than 1. With this threshold, the false alarm rate is 4.0% and the 
missed detection rate is 3.0%. This demonstrates that under the same assumptions, the Bayesian technique 
is extremely effective at detection with quantifiable performance and at shorter integration times. 




D 



log O 10 

Fig. 1. — Value of O™ as a function of Fig. 2. — Histogram of logO" 

n for the two simulated sets, for r = 1. for the two simulated sets, r = 1 
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3.2. Assume I p is uncertain 

The first complication to consider is a lack of knowledge of I p . In Kasdin et al. (2003), we take the 
approach of estimating I p (whether a planet is there or not) and examining the quality of the estimate to 
make the decision. Here, we avoid photometry and still use a Bayesian hypothesis testing approach to make 
the decision in a shorter integration time. 

If we assume under Hi that I p follows a gamma law, T(a, (3), then using Eqs. 3 and 4 we get 

i P " E K- A '''r(a + P ) (py 1 ' a ° 



where the coefficients 7 P (0) are provided in the Appendix. The data have been simulated for r = 1, with 
a = 4 and (3 = Qh/4:. The histogram of the resulting estimated S/N using Eq. 1 and the I p selected from 
the distribution from each run is depicted on Figure 3. Following Eq. 1, if a conclusion is reached only for 
S/N > 5, the missed detection rate reaches 72.7%. As described in Section 3.1, the mean value of I p was 
chosen to correspond to a S/N of 3.52 for r = 1. When I p is sampled from the gamma distribution here 
(with the same mean value), roughly 24% of the time a planet occurs with a mean S/N > 5 that would 
assure a detection by the conventional approach. For the same simulated data, our criterion provides a 
missed detection rate of only 1.0% with the odds ratio decision threshold set at one. 




0123456789 10 

S/N 



Fig. 3. — Histograms of the estimated 
S/N when the planet intensity I p fol- 
lows a Gamma law. The blue boxes 
correspond to the case where S/N < 5, 
i.e. that an existing planet has not been 
detected. 



p(a„|H,) 




p(a,|H ) 

Fig. 4. — ROC curves as a func- 
tion of the relative integration 
time r: the power of the test in- 
creases with t 



The Receiver Operating Curves for this case is depicted in Figure 4. It shows the trade-off between 
the false alarm and the missed detection rates for different values of the relative integration time r. As 
expected, the greater r, the more powerful the test. Nevertheless, note that a relatively high correct decision 
probability p(aa\Ho) + p(ai\Hi) (found via the Monte Carlo results with a decision threshold of one) can 
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be achieved even for low integration times (see Table 1). In practice on real data we may choose the lower 
optimal threshold deduced from the simulations, thus achieving slightly higher probabilities of success. 



T 


0.1 


0.2 


0.5 


1 


p(a \H ) + piailHi) 


0.732 


0.804 


0.876 


0.900 



Table 1 : Evolution of the correct decision probability with r 



3.3. Assume 1^ is also unknown 

The final, and most complicated, case is to assume that both I p and It are unknown. Again, in our 
previous report and in Kasdin et al. (2003), we resorted to two-parameter least-squares estimation (matched 
filtering) on the photometric data. By simultaneously estimating the irradiance of the planet and the 
background, a decision can be made on the planet, but the integration times become quite long. With the 
Bayesian approach, we found that we have a very high probability of success even at the same integration 
time as the previous case. 

For the sake of brevity, we defer the derivation to the Appendix. If p(Ib) = T(a, 6), we use the marginal- 
ization law Eq. 4 to eliminate its dependence from the odds ratio and obtain: 



y^+K 

o» - 1 1 J i>+ P ) (py T{a K - P ) 



where 



A similar Monte-Carlo study as before shows that we can reach a correct decision probability of 0.84 for the 
same integration time ( r = 1) with -§[^j = 1 and E[I%] = E[I b ]/2, E[I%] = E[I b }/2. 



4. Conclusion 

Bayesian hypothesis testing is a promising approach to planet detection and to associate a quantitative 
confidence level with the detection system, with a relative integration time smaller than t\ and without per- 
forming photometry. Moreover this approach allows us to take into account the uncertainty associated with 
the nuisance parameters such as the intensity of the background via the marginalization law. Nevertheless, 
choosing a prior may strongly influence the quality of the result. Further work aims to develop this same 
approach with other priors and to study its robustness. 

No speckle has been considered yet. Current work aims at defining a three- hypothesis test to detect if 
the pixel under study contains a planet, speckle, or just noise. A trio of three bandpass filters, while allowing 
us to consider quasi- monochromatic signals (and directly applying this technique), would also permit us to 



define speckle as a space-variant component. Note that adding other features in addition to just the PSF 
shape in a Bayesian or optimal classifier could also improve the results. 
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A. Appendix 
A.l. PSF Truncation 

Let A be the FWHM of the PSF. The maximum number of pixels that can be covered by the PSF is 
an even number denoted by 2K + 1 where 

As 

p(z k \M?)=p(zk\M ),k$ [n-K,n + K] 

Then 

M n+K 

\[p{z k \M?) n P{*k\M?) 

n n _ k=l = k=n-K g_ £q {lp,lb,0) 

M n+K L n (J. ) 

Y[p{z k \M%) n p{zh\ms) 0Kb > 

fe=l k=n-K 

Example. If K = 1, 5 pixels only can be involved, 

Q n ^ P (Zn-l\M?)p(z n \M?)p(z n+1 \M?) 

10 p(z n ^\M^p(z n \M^p(z n+1 \M^) 
In the following, we will use the following notation 

k=n+K n+K 

o~k — ^ z k ;ir k = ][ z k l 

k —n—K k —n—K 



A. 2. Ip is assumed to be known 

Obviously, as is assumed to be known 

p(Z\M n ) = [ p(Z\M^,I b )p(I b )dI b =p(Z\M^I b ). 
For the constant-rate model Mq 



L (h)= n p(*k\Mz,i b )= n ^ A ? exp (- Afe ) 

k—n — K k=n—K 

n+K 

] f — - (I b AaT) Zk exp (-J 6 Aar) = — (/ 6 tA Q ) m exp (-/ 6 tA5) , 



k=n-K Zkl WK 
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where AS = (2K + 1) Aa is the area of integration around the n'th pixel. 

Assuming that 9 is uniform on [£„, Cn+i], the likelihood of the alternative model M{ 1 is 

rC 



Li(I p ,I b ) = J^ +1 p(Z\M?,I b ,I p ,0) P (0)d9 

1 /-Cn+l «+f 

= Aa / d6 n p(«fc|Mr,/ 6 ,J p ,fl) 

^ n k—n — K 

1 /-Cn+l 1 ™+5 

= a-/ d0 — 1 (4Aar + / p AP fc (0)) 2fc exp(-r(/ b Aa + / p AP fe (0))) 

Compute 

n+K 

M(6) = [] (AAar + rJpAflt^^expC-T^Aa + JpAilb^))) 

k=n-K 

n+K , TAP (ff\ \ Zk ( k=n+K \ 

= (/ 6 rAa)«exp(-J 6 rAS) [] 1+ P , A ) exp ~ r E V^M*) 

fc=n-K V lbT£Xa ' \ k=n-K J 



n-t-i\ , 

= {I b TAaf K exp(-r(/ 6 A5 + / p AP)) [] (l + Q 



"" ' AP k (0) 
Aa 



k=n-K 

k=n+K j 

by noting that AP& (0) = AP and by noting Q = -f-. We then have 

k=n-K b 



U (I P , h) = {I "l^ K cxp (-r (7 6 AS + J P AP)) J " +1 dtf f[ f 1 + Q- 



Aa 



U (h) 



k=n-K 
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A. 3. Assume I p is uncertain 



As h is still assumed to be known, L is unchanged. Using the marginalization law, we have 
p{Ip)Li (Ip)dl p 



o n w (h) 



Lo (4) 



by posing 



For the sake of simplicity, we only present here the case for K = 2 (case described in the main document). 
Using the binomial law 



a+^EffV 



n=0 

2K + 1 times we have that 

n+K , \u fa\\ z k °K 



k=n-K V 7 p=0 v 0/ 



where 



7„w = e 



z n+1 \ /AP„ +1 (6») 



r-max(0,p-2 n ^i) 
min(z n _ K +z n _s:-|_i,r) 



Aa 



p—r 



r (0) = E ai(0,n,n + l)a p _ ( (6»,n + 2,n + 3) 

/— maxfO,r — 2 n -2n + i) 

a <^) = E L z_m Hd 2 \-t~ 



m— mayi(0, 1 — Zj ) 



so that 



1 1 r^+i Z" 00 / r 
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Remark. ./Vote tftaf the computation of j Q Aa J p (0) d6, that may look very time-consuming, can be done 
before any data are collected. To be computed, only the sums (depending on Zi) have to be computed online. 

A. 4. When lb is uncertain 

We now have to apply the marginalization law for both sides. 



/>oc 

Lo = / p(h)L (I b )dI b 
Jo 

\<?K r oo 



r(a)b»n K J '» 6X H Vj 6 

(TAa)' TK T(a')b' a ' 
tt k r (a) o a 
(rAa)" K T (a + o*) o' a ' 



where 



Compute now L\ 
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where 

a" = <jk + a — p = a' — p; 

so 

o" - 1 ( 1 \V^y r( Q +p)rK- P ) 



